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We study the stability of anti-de Sitter (AdS) spacetime to spherically symmetric perturbations 
of a real scalar held in general relativity. Further, we work within the context of the “two time 
framework” (TTF) approximation, which describes the leading nonlinear effects for small amplitude 
perturbations, and is therefore suitable for studying the weakly turbulent instability of AdS— 
including both collapsing and non-collapsing solutions. We have previously identihed a class of 
quasi-periodic (QP) solutions to the TTF equations, and in this work we analyze their stability. 
We show that there exist several families of QP solutions that are stable to linear order, and we 
argue that these solutions represent islands of stability in TTF. We extract the eigenmodes of 
small oscillations about QP solutions, and we use them to predict approximate recurrence times 
for generic non-collapsing initial data in the full (non-TTF) system. Alternatively, when sufficient 
energy is driven to high-frequency modes, as occurs for initial data far from a QP solution, the TTF 
description breaks down as an approximation to the full system. Depending on the higher order 
dynamics of the full system, this often signals an imminent collapse to a black hole. 


I. INTRODUCTION 

Of the maximally symmetric solutions of the Einstein 
equation, nonlinear stability in general relativity has 
been proven for both Minkowski [T] and de Sitter [5] 
spacetimes. In contrast, the question of stability of anti- 
de Sitter (AdS) remains formally open. A key differen¬ 
tiator of AdS as compared to its A > 0 counterparts is 
that with non-dissipating boundary conditions at infinity, 
perturbations cannot decay and energy is conserved [3]. 
Based on knowledge of nonlinear wave propagation in 
the absence of dissipation, AdS has been conjectured to 
be unstable [H [5] (see also 0)- This expectation has 
been corroborated by numerical simulations—supported 
by perturbative arguments—which showed that certain 
initial configurations evolve to black holes, no matter how 
small the initial deviation from AdS was taken [7]. This 
work showed, further, that the eventual gravitational col¬ 
lapse resulted from a turbulent cascade of energy to high- 
frequency modes of AdS, mediated by resonant interac¬ 
tions. 

Instability of AdS would have implications for a num¬ 
ber of fields, ranging from potential gravitational insta¬ 
bilities in other low-dissipation or confining geometries, 
to thermalization of conformal field theories (CFTs). In 
the context of AdS/CFT (within the regime where gen¬ 
eral relativity holds in the bulk), the formation and sub¬ 
sequent evaporation of a black hole in AdS is believed 
to be dual to the process of CFT thermalization. The 
more recent discovery of initial configurations in the bulk 
that appear to avoid black hole formation [5] was, there- 


* sgreen@perimeterinstitute.ca 
^ antoine.maillard@ens.fr 
^ llehner@perimeterinstitute.ca 
§ steve.liebling@liu.edu 


fore, somewhat surprising, as that would indicate non- 
thermalizing CFT configurations. This finding led to the 
identification of several “islands of stability” in AdS 0- 

m- 

Significant progress towards an analytic understanding 
of the dynamics was achieved with the introduction of a 
powerful perturbative framework—the two time frame¬ 
work (TTF)—for analyzing small perturbations of AdS in 
terms of coupled nonlinear oscillators [12] (see also [I3j). 
This framework efficiently captures the resonant energy- 
exchange interactions between normal modes, while effec¬ 
tively “integrating out” high-frequency oscillations. TTF 
led to the discovery of a pair of quantities—the energy E 
and particle number N —that are conserved at the lead¬ 
ing nonlinear level m nnj. These quantities play key 
roles in understanding long-term dynamical behaviors, 
including dual (direct and inverse) turbulent cascades 
and non-equipartition of energy m- 

The main purpose of this paper is to establish a 
large new class of islands of stability within the TTF 
approximation. The central stable equilibria—quasi- 
periodic (QP) solutions [m [TS] —form discrete families, 
each family itself parametrized by N and E. In this 
work we (i) construct the families of equilibrium so¬ 
lutions, (ii) perform a linear stability analysis within 
TTF showing stability, and (iii) use the results to un¬ 
derstand the long-term behavior of both collapsing and 
non-collapsing initial configurations. In particular, the 
stability analysis gives rise to a perturbation spectrum 
that agrees with and explains “recurrences”—long-term 
nearly-periodic approaches of the configuration to the ini¬ 
tial state, first observed in the Fermi-Pasta-Ulam (FPU) 
system of coupled oscillators [TS|— which were observed 
numerically in the full system. Dependence of the fami¬ 
lies of QP solutions on the two continuous parameters N 
and E extends previously known one-parameter families 
(time-periodic solutions (TU]) and provides a clear con- 



2 


nection between conserved quantities and stable islands. 

A. Background 


At lowest nonlinear order, we showed that gravitational 
self-interactions of (j) are taken into account provided the 
coefficients Aj{T) satisfy the coupled ordinary differential 
equations m, 


Following [7], we restrict analysis to the spherically 
symmetric case and four spacetime dimensions. As a 
proxy for gravitational degrees of freedom, we take as our 
model a real massless scalar field (j) coupled to general rel¬ 
ativity. Ignoring gravity, the scalar field is characterized 
by normal modes with spatial wavefunctions 


ejix) = 4 


+ cos^ X 3 + J; I sin^ 

( 1 ) 

and frequencies ujj = 2j + 3 (j = 0,1, 2,...). 

Since the frequency spectrum is commensurate, non¬ 
linear gravitational interactions are resonant, and those 
interactions cause energy to be readily transferred among 
the mode^ Numerical simulations have shown that for 
certain initial data, energy is transferred from low-j to 
high-j modes—a direct turbulent cascade [7]. This cas¬ 
cade concentrates the energy—the high-j modes are more 
highly peaked in position space—and eventually leads to 
black hole formation. The cascade behavior persists self- 
similarly as the amplitude e of the initial scalar field is 
decreased, with the time to collapse scaling as 1/e^ (see, 
e.g.. Fig. 2 of H). 

In contrast, other initial data seem to avoid collaps^ 
as e —>■ 0. In addition to direct cascades, these solutions 
feature inverse turbulent cascades, which transfer energy 
to low-j modes [HHI]. Collapse is avoided if the inverse 
cascades sufficiently hinder the flow of energy to high-j 
modes. 

A key observation is that energy cascades and nor¬ 
mal mode oscillations are governed by independent time 
scales. As e —^ 0, nonlinear interactions become weaker— 
the stress-energy tensor xt oc e^, and gravitational self¬ 
interactions of j) scale as —so the energy transfer time 
scale is proportional to 1/e^. Meanwhile, normal mode 
oscillations proceed independently of e. This separation 
of time scales means we can use multiscale analysis meth¬ 
ods to study the slow mode-mode interactions indepen¬ 
dently of the fast normal mode oscillations in the limit 

e ^ 0 [ig. 

We define the “slow time” r = tje^. Over short time 
scales the scalar field is well-approximated as a sum over 
normal modes. Thus, we take as ansatz cj) = with 

OO 

(Aj(T)e-*“’’‘ -h ej{x). (2) 

i=o 


^ The frequency spectrum is also resonant with a massive scalar, 
for other spacetime dimensions, and in the absence of spherical 
symmetry. 

^ Simulations are of finite duration, and the limit e ^ 0 cannot be 
obtained numerically, so collapse-avoidance is a conjecture. 


(3) 

klm 

known as the two time framework (TTF) equations. The 
TTF equations were also derived using renormalization 
group perturbation methods to re-sum secularly grow¬ 
ing terms that arise in ordinary perturbation theory [13] . 
Notice that the TTF equations possess the same scaling 
symmetry, A{t) —>■ eA(Tfe^) seen in the full (non-TTF) 
system in the limit e —>■ 0. 

The numerical coefficients appearing in (§ arise 
from overlap integrals involving the ej(x), and they van¬ 
ish unless j + k = I + m. This fact, together with the 
specific form of the equations ([^ (i.e., the lack of terms 
such as A^AiAm, etc.), arises because the only resonances 
that are present in the system are those such that m 

U)j -\- Ulk = tVi -\- U)m ■ ( 4 ) 

This property is related to a hidden symmetry in 
AdS [H]. For further discussion on the absence of certain 
resonance channels, see m- 

Within their regime of validity, the TTF equations 
yield approximate solutions much more economically 
than full numerical relativity simulations |12j . Indeed, a 
significant speedup is gained by not modeling the rapid 
normal-mode oscillations. Moreover, the TTF approxi¬ 
mation improves as e —>■ 0—a limit that is especially hard 
to reach in numerical relativity. Nevertheless, in the same 
way that finite difference methods employ a discrete spa¬ 
tial grid, the set of TTF equations must in practice be 
truncated at finite j = jmax (similar to pseudo-spectral 
methods). Previously, we computed (by performing ex¬ 
plicit integrations on a mode-by-mode basis) the coeffi¬ 
cients up to jmax = 47 [ig. We now have closed 
form expressions for the coefficients (see App. 0 that 
enable us to work to much larger jmax- We typically set 
jmax = 200 in this paper, which in many cases provides 
an excellent approximation. In particular, the recurrence 
dynamics of non-collapsing solutions are well-captured. 

While useful as a calculational tool, the main power 
of TTF is analytic]^ Indeed, in m US] it was uncov¬ 
ered that the TTF equations conserve a total of three 
quantities: The total energy and particle number, 


j 

(5) 


(6) 
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^ See also Eq] for another recent illustration of the power of this 
approach within general relativity. 
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as well as the Hamiltoniarj^ 

"i E -^Y.c,\A,\\ (7) 

jklm j 

where Cj are additional constants. Conservation laws of 
E and N are associated with two C/(l) symmetries, 

A,{t)^A,{t)E-^^, ( 8 ) 

A,{T)^A,{r)E<^, (9) 

respectively, for d G M constant; conservation of H is 
associated with time-translation symmetry |14] . (The 
symmetries and associated conservation laws were first 
uncovered for the TTF equations that describe a non¬ 
gravitating scalar field in AdS 4 , with quartic self¬ 
interaction V{(j)) = X(jA/A\ [H].) Simultaneous conserva¬ 
tion of E and N implies that direct and inverse turbulent 
cascades must occur together, and that energy equipar- 
tition is in general not possible m- 

Finally, we showed in m that the TTF equations give 
rise to equilibrium solutions, which are QP. That is, each 
mode amplitude, 

= (10) 

with f3j G M. Simulations in TTF and full numerical rel¬ 
ativity both provided evidence for stability of these QP 
solutions. The case was then made in [TS] that general 
non-collapsing solutions can be treated as perturbations 
about associated QP solutions—in other words, QP solu¬ 
tions with the same E and N. As an example application, 
we studied two-mode initial data, which exhibits FPU- 
like [TB] recurrences over long time scales. We showed, by 
interpolating initial data between two-mode and associ¬ 
ated QP, that the recurrence times were only marginally 
affected. We therefore concluded that a proper stabil¬ 
ity analysis might predict these times, and QP solutions 
might provide anchor points for the “islands of stability” 
in AdS. 


B. Summary 

In this paper we present a comprehensive analysis of 
QP solutions and their relation to AdS (in)stability. Af¬ 
ter presenting the algebraic equations governing QP so¬ 
lutions in Sec.|ITj we show that they extremize E[ for first 


As described in detail in the system § in the “origin-time” 
spacetime gauge of [3[T2l[T3] is not a Hamiltonian system itself. 
However, in the “boundary-time” gauge of [3 [H] the system 
is Hamiltonian with Hamiltonian H. Both gauges possess the 
same conserved quantities, so in this paper we shall refer to H 
as the “Hamiltonian”, despite working in origin-time gauge (for 
comparison with prior numerical simulations). Note also that 
the equivalent expression in m did not include the second term 
in H. 


order variations holding E and N fixed. We then numeri¬ 
cally map out the space of solutions to the QP equations. 
This space can be divided into a number of families of 
solutions, each one depending on two continuous param¬ 
eters, E and N. Because of the scaling symmetry of the 
TTF equations, these families are scale-invariant, so it is 
often useful to exchange E and N for an overall scale, 
and the ratio T = E/N —which we identify with the 
“temperature.” 

In Sec. m we perform a linear stability analysis of 
QP solutions within TTF. We uncover two 2-dimensional 
subspaces of special perturbations: the first corresponds 
to a pair of generators of the C/(l) symmetries ^ and 
of TTF; the second represents infinitesimal perturbations 
to nearby QP solutions with different E and N. (We 
make use of these special perturbations to generate the 
continuous families of QP solutions parametrized by E 
and N in Sec. |llj) The remaining perturbations preserve 
E and iV, and may be decomposed into eigenmodes de¬ 
scribing small oscillations. The corresponding eigenval¬ 
ues determine stability. We present a numerical method 
to perform this stability analysis given any particular 
background QP solution. 

After presenting the framework for analyzing stability, 
we apply it to the families of QP solutions identified in 
Sec.|^ We argue, by explicitly checking a large number of 
QP solutions, that the “physical” families—those that do 
not depend strongly on the mode cutoff jmax in the limit 
Jmax —>■ oo —are all stable. By contrast, QP solutions 
that are not members of these families can have unstable 
modes. 

In Sec |Iy] we apply the results of the stability analy¬ 
sis to understand long-term evolutions. Since E and N 
are conserved, motion in phase space is constrained to 
constant-(A, iV) hypersurfaces. Each surface intersects a 
given stable QP family at most once, resulting in a dis¬ 
crete collection of QP solutions. If initial data lies within 
the i7-trough around one of these QP solutions with the 
same {E, N) , then we associate it to that QP solution. 
Under evolution the solution is then confined to oscillate 
about its associated QP solution. We illustrate, through 
several examples, how nonlinear evolutions of initial data 
within TTF inherit many of the properties uncovered by 
the linear stability analysis. In particular, the linear sta¬ 
bility analysis explains nonlinear recurrences as oscilla¬ 
tions about QP solutions, and the eigenmodes (and com¬ 
binations thereof) predict the recurrence times. Thus, we 
obtain approximate recurrence times without performing 
any time integrations. This approach to understanding 
recurrences is a generalization (to two conserved quanti¬ 
ties) of the g-breather approach to understanding FPU 
recurrences 

For evolutions that remain close to stable QP solu¬ 
tions, the energy spectra remain close to the exponential 
energy spectra of the QP solutions. In contrast, solu¬ 
tions that are not close to QP solutions tend to approach 
power laws in TTF, consistent with earlier studies [23]. A 
power law spectrum contains far more energy at high-j. 
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and when translated to a description involving spacetime 
fields at finite e, the energy is far more concentrated at 
the origin of AdS; in fact, it fails to even converge in 
j. When deviations from AdS become large, TTF no 
longer applies, and higher-order dynamics take over. It 
is often the case that the higher order dynamics rapidly 
drive collapse once they take hold [7]; the role of TTF is 
to indicate whether this regime is reached. 

It is important to keep track of the various levels of 
levels of approximation used in this work, so we summa¬ 
rize them here. First, the TTF equations are taken as 
an approximation to the full system, valid in the limil[^ 
e —>■ 0. Secondly, we truncate the TTF system at a fi¬ 
nite number jmax of modes. Finally, we perform a linear 
stability analysis of QP solutions within the truncated 
TTF system. Throughout this work we will address the 
validity of the various approximations. 

II. QUASI-PERIODIC SOLUTIONS 

There is already strong evidence that there are sta¬ 
ble equilibrium solutions— islands of stability —in AdS, 
namely the time-periodic solutions [3 [TO]. These solu¬ 
tions are nonlinear generalizations of individual normal 
modes, with the effect of gravity being to shift the fre¬ 
quency. Such solutions are moreover realized as solutions 
to the TTF system ([^ of the form 

A,(t) = (11) 

for some fixed mode number k. (The analysis of [TO] , 
however, is accurate to higher order in e.) For a 
given k there exists a 1-parameter family of solutions, 
parametrized by Afc(O), or equivalently, the energy E. 

Inspired by the periodic solutions, we identified in |12j 
a much larger class of gitasz-periodic (QP) solutions. Al¬ 
lowing for all modes to be excited periodically (but with 
different periods), we sought solutions of the form 

Mr) = 

with {aj,fij) G C X K. Such solutions would have con¬ 
stant energy, Ej, in each mode—finely tuned so that en¬ 
ergy flows between modes are perfectly balanced. 

Substituting the ansatz above into ([^, we have 

klm 

( 12 ) 

We see that the r-dependence may be canceled from both 
sides by imposing the condition 

Pj = Po + {Pi - Po)j, (13) 


^ For an interesting discussion of when solutions of the approxi¬ 
mated system might correspond to solutions of the full system, 
see |25| . 


reducing the system to 

- 2ujj [/3o + j (/3i - /3o)] aj = ^ sl^l^ataiam- (14) 

klm 

Without loss of generality, henceforth we take aj G M 
in the equation above (this represents a choice of ini¬ 
tial time r = 0). We thus have jmax + 1 algebraic 
equations for jniax + 3 unknowns. That is, we have 
two free parameters—one more than the time-periodic 
solutions—which we will often take as E and N. 

A. Extremization of H 

Quasi-periodic solutions extremize the Hamiltonian E[ 
with respect to perturbations that preserve E and N . To 
see this, we first introduce some notation (following [2]). 
We split the coefficients 

‘^kim ~ ‘^jklm + '^jk {^jlPkm + Sjmdkl) , (15) 

where 5®^;^ is symmetric under interchange of jk with 
Im (as well as exchange of j with k or I with m). The 
quantity is antisymmetric and takes the form 

n%=cMk-Cku:l (16) 

We then define the quantity 

V = sjfil^AjAkAiAm- (17) 

jklm 

It may be shown that 

gy _^ 

~ 2 ^ { ^iklmAkAlAjn- (18) 

^ klm 

Thus (§ may be re-written 

-2iujj-^ =-^+ ‘^'^'^fk\Ak\‘^Aj, ( 19 ) 

^ k 

or in terms of the Hamiltonian, 

iujj-^ = -I- 2uj'jAj (20) 

^ k 

Note that the presence of the last term indicates that the 
system is not actually Hamiltonian in the “origin-time” 
spacetime gauge in which we work (see footnote . This 
term is not present in the “boundary-time” gauge [14] . 
Now consider a variation that fixes E and N, 

=E + M’) 

j ^ 

- - {AjSAj + AjSAj) 

jk 

(^“j - 

j ^ '' k 
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On the second line we used the TTF equation (20). On 
the last line the final term vanishes for variations that 
preserve E. For Aj also quasi-periodic, we can now use 
the ansatz ( [To| ) and to simplify the first term, 


SH — ujjPj {AjSAj + AjSAj) 

3 

— X] ^3 [/^o + 3{Pi - Po)] {AjSAj + AjSAj) 

3 


= 0 , 


3o - Po) 


SN 


( 22 ) 



since we fix E and N. Thus, QP solutions are critical piG. 1. Energy spectra for several QP solutions that were 
points of E[ for perturbations that fix E and N. obtained numerically. These solutions are all members of the 

jr = 0 family. Here we take jmax = 40. 


B. Families of solutions 


The QP equations (14) have two free parameters, 
which must be fixed prior to solving. But, even after do¬ 
ing so, there remain multiple solutions because the equa¬ 
tions are nonlinear. This gives rise to multiple families of 
QP solutions, each extending over some range of E and 
N. 


We solve the QP equations numerically, following sev¬ 
eral approaches described in App.|^ As always, the TTF 
system is truncated at j = jmax < oo, and the physical 
continuum limit corresponds to jmax —>■ oo. Thus, any 
QP solution that depends strongly on jmax must be dis¬ 
carded as unphysical. 

The simplest way to obtain QP solutions (used in [T^ l 
is to use a Newton-Raphson method, which works well 
if a good initial seed can be chosen. Since we know that 
single-mode configurations (111 are solutions, we search 
for solutions dominated by single modes j = j^, but that 
have nonzero energy in the other modes. The energy 
spectra Ej = of several such solutions from 

the jr = 0 family are illustrated in Fig. Rather than 
parametrizing the solutions by the continuous parameters 
E and N we have labeled the spectra by the temperature 
T = E/N. The other parameter is simply an overall scale 
that does not affect the shape of the curves. 

Notice that for small T, the energy spectra approach 
exponentials. (The minimum temperature for the jr = 0 
family occurs in the single-mode limit, with Tmin° = 


Eq/Nq = ujq = 3.) For larger T the spectra deform 
and it becomes increasingly difficult to obtain solutions 
using the Newton-Raphson method. For such cases we 
can obtain solutions by perturbing known solutions to 
different E and N (see App. E- For the jr = 0 fam¬ 
ily, solutions exist up to T = T^ax = = 2jiaax + 3, 

which is the maximum possible temperature for the trun¬ 
cated collection of modes. Such solutions are highly 
deformed from exponentials—the maximal solution has 
all energy in mode j = jmax—and are not physical be¬ 
cause of the dependence on mode truncation. Requiring 



FIG. 2. Energy spectra of QP solutions from several discrete 
families (with different jr). The temperature in each case is 
very close to as the QP solutions shown here are very 
close to single-mode solutions, (jmax = 100) 


T Tmax will select for physical configurations, and, 
with this restriction, the physically relevant spectra are 
all nearly exponential. Extrapolating to the continuum 
limit jmax oo —where by definition there are no un¬ 
physical solutions—we expect all jr = 0 solutions to have 
nearly-exponential spectra (for any T). 

In Fig. we plot the spectra of QP solutions from 
families with various jr > 0. Each solution is peaked 
s-t j = jr, and decays exponentially to both sides (with 
slight deformation for j < j,.). As jr increases, so does 
the minimum temperature Tmin ~ respective 

QP family. We find that the jr > 0 families do not 
extend in temperature all the way to Tmax (in contrast 
to the jr = 0 case), but that the range of temperatures 
increases with jmax (see Fig. [^. In the jmax —>■ oo limit 
it is not clear whether the families have a finite or infinite 
extent. 

It is possible to construct additional QP families. For 
example, the resonance condition @ implies that if 
only even-numbered modes are excited initially, they will 
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FIG. 3. The domain of existence of QP families for jV € 
{0,1, 2,3, 4, 5} and jmax £ {10, 30, 50}. For jV = 0, the family 
is defined in the full domain [3,2}max + 3]. Note that the 
bounds of the vertical axis increases with jmax. 



FIG. 4. QP solution not smoothly connected to single-mode 
solution, (jmax = 100) 


never excite odd-numbered modes. In this case QP so¬ 
lutions can be found that are similar to those of Fig. 
but skipping every other mode. Finally, there are so¬ 
lutions that have considerable energy in high-j modes 
that do not appear to connect to the families above (see 
Fig.i. These latter solutions are clearly dependent on 
mode-truncation, so we discard them as unphysical. 


III. STABILITY OF QUASI-PERIODIC 
SOLUTIONS 

In [12], we numerically tested the stability of several 
QP solutions in the jV = 0 family within the full (non- 
TTF) theory. For the duration of the simulations, per¬ 
turbations oscillated about the QP solutions over time 
scales long compared to the AdS crossing time. 

To address stability more systematically, in this sec¬ 
tion we undertake a linear stability analysis of QP so¬ 
lutions within TTF. In Sec. |III Aj we linearize the equa¬ 
tions (|^ about an arbitrary background QP solution. We 


show that through an appropriate change of variables, 
the time-dependence (resulting from a time-dependent 
background solution) can be eliminated, leaving an au¬ 
tonomous system of the form. 


dx 

dr 


= Ax. 


The matrix A depends on the background QP solution, 
and is independent of time r. The problem of solving in 
time for the perturbation vector x is, therefore, equiva¬ 
lent to that of diagonalizing A. 

In Sec. IIIB we identify special solutions (infinitesi¬ 
mal U{1) symmetry transformations and perturbations 
to other QP solutions) unrelated to stability, and in 
Sec. |III C| we outline the numerical procedure for find¬ 
ing the remaining eigenvalues of A. Finally, in Sec. HID 
we apply this approach to study the stability of the fam¬ 
ilies of QP solutions identified in Sec. m We sample 
a large number of solutions within the “physical” fami¬ 
lies, and find that they are all Lyapunov stable—initially 
small perturbations remain small, but they do not decay 
to zero (see Sec. 23 of (Hj). In Sec. HIE we comment on 
nonlinear stability for finite-sized perturbations. 

Throughout this paper, we work in the origin-time 
spacetime gauge (see footnote]^. We note that the sta¬ 
bility analysis would go through nearly identically in the 
boundary-time gauge, where the system is truly Hamil¬ 
tonian. The change of gauge only contributes a time- 
dependent phase shift to the TTF coefficients A^ (r), and 
our stability results hold in both gauges. 


A. Linearized equations 


Consider a perturbed QP solution, 

+Cj(T), 


(23) 


where A^^(t) = a^-e with {aj,f3j} G K. Substi¬ 
tuting into ^ and keeping terms to first order in fj, we 
have 




dr 2uj 


klr 


(24) 




Since the background QP solution has quasi-periodic 
time-dependence, so do the coefficients of this equation. 
If the coefficients were in fact periodic one could have ap¬ 
plied the Floquet theory to obtain the general solution 
to (24) in terms of eigenmodes, and thereby determine 
stability (see Sec. 28 of [26]). (In fact, by tweaking the 
values of /3o and /3i so that they are rational multiples 
of each other, periodicity can be achieved, although the 
period might be quite long.) The Floquet approach re¬ 
quires numerical integrations over one period to identify 
the eigenmodes, which is somewhat tedious, but works 
generically for periodic systems. 
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For our TTF system, however, the analysis simplifies 
due to the resonance condition. First, factor out the 
background time-dependence in each perturbative mode 
to define new variables, 


This gives rise to the autonomous equations 

dT 


(25) 


(26) 


respectively. 

It is straightforward to check that these perturbations 
satisfy (28)-(29). Indeed, (29) holds trivially, while (28) 
holds because of the resonance condition 0 [in the case 
of (33)] and the QP equation (14). 


+ ^ X! (Xkaiam + akXiam + akUix-m) ■ 


klm 


These equations contain complex conjugations of Xj 
are therefore not linear over C. To obtain a linear system, 
split Xj into its real and imaginary parts, 

Xli^) =Uj{T)+iVj{T). (27) 

The system is now reduced to 
dui 


2. Perturbations to nearby QP solutions 

Consider now a perturbation from a QP solution to 
another QP solution, 

^ (a^. + ( 35 ) 

The new QP solution is required to satisfy the QP equa¬ 
tion (14) as well. To first order in the perturbation, this 


requirement takes the form 
—2u)j {ttjSBj + Bj^oij) 
:U) 


(36) 


dr 


= -Pj^j 


(28) 


^ ^ Bklm OLkOti6(y.jYi) j 


klr, 


2uj 


^ ^ ^ ^klm ( ^l^m'^k ^m^k'^l “ 1 “ ^l^k'^m) 7 


klm 


and the condition (13) implies either of 

S^j —y (jjjO^ 


^'^0 a 
dr 

1 


(29) 


SBj —>■ 0, 


(37) 

(38) 


for 9 e 


+ E infinitesimal version of the perturbation is 

Zcj - ... 


klm 


the linearized energy, particle number, and Hamiltonian, 


It can be shown that the equations (|28|)-(|29|) conserve 

nian, 

(30) 

(31) 




Saj 
-ajrSBj 


(39) 


SE = 8Y, 


LOjOjUj, 


Using this mapping it is easily checked that (36) is iden¬ 
tical to ([ 2 ^, and that (1^ holds for both cases ([3^ 


and (38). 


SN = 8'^ujjajUj, 


6H=- /?i-/3o-4^C,ad JU 


--(5/3o-3/3i)5fV. (32) 


B. Special solutions 

1. U{1) symmetry transformations 

Recall that the TTF equations are invariant under two 
U{1) symmetries 0-0, 

A,{t) ^ A,(r)e-^^ 

for 0 G M constant. Off of QP solutions, infinitesimal 
t/(l) transformations take the form 


Perturbations (371 and (38) represent a 2-parameter 


family of solutions to the linearized equations. This fam¬ 
ily can be re-parametrized in terms of SE and 6N, allow¬ 
ing for the families of QP solutions in Sec. |IIB| to be fully 
obtained as orbits of these perturbations (see App. H- 

Together, infinitesimal 17(1) transformations and in¬ 
finitesimal perturbations to nearby QP solutions form 
two 2-dimensional generalized eigenspaces (with eigen¬ 
value 0) of the matrix A representing the linear system 
(see below). Indeed, the action of A on a perturbation 
of the form (37) gives a 17(1) transformation (33), and a 


(34) 


subsequent action of A gives 0. [Similarly, (38)- 
0.] A 2-dimensional generalized eigenspace does give rise 


to linear growth in the solution [see (39)], but this growth 
is not relevant to the question of stability since it is sim¬ 
ply an infinitesimal perturbation to another equilibrium 
solution (35). 


C. General solution technique 


Uj 

Vi 


Vi 


0 

ujjaj9 


aj0 


(33) 

(34) 


It is convenient to express (28)-(29) in matrix form. 
Defining 


X = 


(Uj) 

W) 


(40) 




































the perturbative equations take the form 


dx 


( 41 ) 


where A is a (2ji„ax + 2) x (2jniax + 2) constant real 
matrix. We now complexify the equation and put A in 
Jordan form, taking real solutions in the end. 

In general, the background QP solution, and hence the 
matrix A, are known only numerically. This is problem¬ 
atic since the Jordan decomposition is numerically ill- 
conditioned—if A has multiple eigenvalues, small errors 
in A can lead to large errors in its Jordan form. In par¬ 
ticular, we know from the previous subsection that A has 
two generalized eigenspaces of dimension 2, which can be 
misidentified as distinct 1-dimensional eigenspaces. 

In contrast to the Jordan decomposition, the Schur 
decomposition is well-conditioned numerically and con¬ 
tinuous in the matrix elements. We therefore perform a 
Schur decomposition of A, 

S = U-^AU. (42) 


Here U is unitary, and the matrix S is upper triangu¬ 
lar with eigenvalues along its diagonal. The known gen¬ 
eralized eigenspaces of A have eigenvalue 0. Numeri¬ 
cally, however, these may deviate slightly from zero. We 
also find that, generically, all other eigenvalues are well- 
separated from 0. So, to correct the errors in the gener¬ 
alized eigenspaces, we round off all infinitesimal diagonal 
components of S to 0, and denote this new matrix S. 

Finally, we take the Jordan decomposition of S, 


, -yi -^-1-^-1-1— 

- ^ - 

Alo - 

:r ■ 

As 

. 


-a- Ae <* 
A5 >< 
0 A4 

:r )’ - 

-a- A2 , 
-a- Ai 

_^_ 

°0000000000000000000000000000 

.. 

*000000000000000000000000000000 



_1_^_1_^_1_^_ 


15 20 25 30 35 40 45 50 


Jmax 


FIG. 5. Dependence of eigenfrequencies on truncation jmax. 
We plot the ten lowest eigenfrequencies for the QP solution 
with T = 3.75 and E = 8. After decreasing noticeably up to 
jmax ~ 25, eigenfrequencies approach asymptotic values. 


While our system is not Hamiltonian (in the “origin¬ 
time” gauge used here El), the pattern of eigenvalues 
is nevertheless preserved. In particular, for any decaying 
mode there exists a corresponding growing mode, so the 
best one can hope to achieve in terms of stability is Lya¬ 
punov stability. In this case, all modes have harmonic 
time-dependence with no growth or decay—i.e., purely 
imaginary A. 


D. Results 


J = p-^'sp. 


(43) 


The matrix J always contains the expected pair of 2 x 2 
Jordan blocks. Aside from these, we found that the Jor¬ 
dan form J was always diagonal. We denote the addi¬ 
tional (2jmax — 2) eigenvalues by A„, and the associated 
eigenvectors of A (the column vectors of UP) by e„. 

To obtain the time-evolution of a linearized perturba¬ 
tion of a QP background (of the same E and N) one must 
project initial data onto the eigenvectors {e„}. Each of 
these eigenvectors then evolves independently as 
If the initial data is real then a real solution is guaran¬ 
teed. 

There are relationships between the eigenvalues of A. 
Since A is real, if A is an eigenvalue, then so must be 


A. Also, since A is of the form 


0 -C 
D 0 


[see (281- 


(29)], A^ = 


( CD ) ’ eigenvalue of A^ 


occurs twice. Since these eigenvalues are the squares of 
eigenvalues of A, and the eigenvalues of A (excepting 0) 
are generically non-degenerate, if A is an eigenvalue of A, 
then so must be —A. In sum, (A, —A, A, —A) must all be 
eigenvalues. 

These properties of the eigenvalues are also charac¬ 
teristic of symplectic flows and a Hamiltonian structure. 


We applied the above analysis to a sampling of the QP 
solutions described in Sec. nm In almost all cases we 
found that all of the eigenvalues A„ were purely imagi¬ 
nary, implying stability. The only unstable QP solutions 
were those previously deemed “unphysical”, as in Fig. [^ 
and in these cases only a small number of eigenvalues 
had nonzero real part. Therefore, we expect that all of 
the “physical” QP solutions are stable. For these stable 
solutions, we denote the conjugate eigenvalues by using 
negative indices, A_„ = —A„. 

We studied the dependence of the eigenvalues on jmax- 
As jmax is increased by 1, a pair of higher frequency (con¬ 
jugate) eigenmodes is introduced, while (the norms of 
the) existing eigenvalues are shifted slightly lower. In the 
continuum limit jmax —>■ oo, the eigenvalues appear to ap¬ 
proach asymptotic values (see Fig.[^. In that sense, the 
behavior of the low-frequency modes is robust to mode- 
truncation. 

Of particular interest is whether the frequency spec¬ 
trum is itself resonant, as this may imply chaotic dy¬ 
namics at the nonlinear level. In fact, at high frequen¬ 
cies the separation between subsequent eigenmodes A„ 
approaches a constant value as 


iXn — Cl C2n O 



(44) 
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|An| 


FIG. 6. For the QP solution with T = 3.75 and E = 8, we 
plot the magnitude of the components (e„)j of the linearized 
eigenvectors, as a function of eigenfrequency iXn- We see that 
low-frequency eigenvectors trigger low-j modes, (jmax = 50) 



where Ci and C 2 are constants depending on the par¬ 
ticular QP solutioij^ Thus, the high-frequency part of 
the spectrum approaches a commensurate spectrum only 
asymptotically in n. (We will see later that C 2 is closely 
related to the recurrence time for non-collapsing solu¬ 
tions.) 

For perturbations of QP solutions, it is also instruc¬ 
tive to examine the overlap between the original nor¬ 
mal modes of AdS (j-modes) and the QP eigenmodes 
(u-modes). The generic solution to (41) is 




CjiC 


(46) 


where the c„ are constants. For each j, Fig. plots 
the components (e„)^- as a function of eigenfrequency 
iXn- This shows that for initial perturbations consist¬ 
ing of low-j modes, low frequency n-modes are excited. 
Conversely, low-n eigenmodes excite low-j normal modes 
most strongly. This observation explains why low-j 
modes are typically seen to oscillate with the lowest fre¬ 
quencies (see, e.g.. Figs. 10 and 12). 


E. Nonlinear stability 

The linearized analysis above provides useful informa¬ 
tion and intuition for finite-sized deviations from QP so¬ 
lutions as well. Since E and N are conserved quantities, 


For single-mode solutions with jr = k, the eigenvalues may 
be computed analytically, 

*A„ = [A,(0)]^ (45) 

\LJk CJn J 

provided Afc(O) E M (consistent with previous results showing 
stability [El). From this spectrum, the expansion \iA\ may be 
checked explicitly, and the constants Ci and C 2 computed. 


FIG. 7. Hamiltonian plotted as a function of two parameters, 
Hi and fj, 2 , that interpolate [within a constant-(i?, N) surface] 
between 3 QP solutions with T = 7.1. The three QP solu¬ 
tions (blue dots) are members of the jr = 0,1, 2 families. The 
normalization Hq is the value of H for the jr = 0 QP solu¬ 
tion. Note that the ridges result from choosing a non-smooth 
interpolation. 


motion in phase space is constrained to constant-(iil, N) 
hypersurfaces. Each of these surfaces, in turn, intersects 
the families of stable QP solutions at most once each. 
Thus, given the temperature T of the initial data, there 
is a finite number of potentially-relevant QP solutions, 
and these can be determined from Fig. 

Within a given {E, iV)-surface, we know that the QP 
solutions extremize the Hamiltonian H, which is also con¬ 
served in time. In App.[C]we show that, in fact, the stable 
QP solutions minimize El. Since H is conserved in time, 
the size of the surrounding valley in H determines the size 
of the island of stability of the QP solution. One could in 
principle check to see whether given initial data lie within 
one of the valleys, in which case they would remain near 
the QP solution indefinitely (within TTF). As we will see 
in the following section, nonlinear solutions often depend 
closely on the properties (such as the spectrum {A„}) of 
linearized perturbations about QP solutions. 

It is instructive to visualize the minima of H. In 
Fig.0 we plot the value of H over a two-dimensional 
slice of a constant-(£', N) surface. The slice was cho¬ 
sen to pass through three minima, corresponding to sta¬ 
ble QP solutions. Note, however, that the full prob¬ 
lem has a large number of dimensions, with (2jiiiax)“ 
real-dimensional constant-(i?, N) hypersurfaces. More¬ 
over, the continuum limit takes jmax —t 00 . While a 
valley within a finite-dimensional space must have a fi¬ 
nite size, it is possible that this size asymptotes to zero 

^.S Jmax ^ 
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IV. APPLICATION TO ADS (IN)STABILITY 


We now return to our main questions: for a self- 
gravitating scalar field in AdS, how can we predict which 
initial data will collapse in the limit of e —)■ 0? How do 
recurrences arise? How does collapse in the full Einstein- 
scalar system connect to behavior in TTF? 

The TTF equations provide a good approximation if 
the amplitude of the AdS perturbation is small so that 
normal-mode oscillation time scales and mode-mode en¬ 
ergy transfer time scales decouple. The approximation, 
therefore, always breaks down prior to black hole forma¬ 
tion. Knowledge of this fact alone, however, indicates 
that a great deal of energy has transferred to high-j 
modes, and in many cases, subsequent evolution will lead 
to collaps^ 

For small, but finite, perturbations, a simple criterion 
for checking whether TTF has broken down is to evalu¬ 
ate spacetime quantities {4>,gab) from {Aj^r)} and check 
for black holes (i.e., check whether the metric quantity 
A of [7] vanishes at any point, or whether the energy 
in the scalar field satisfies the “hoop-conjecture” [29]). 
Likewise, one could check whether {dt(j>)'^ becomes large. 
We will refer to the blow-up of spacetime fields as “col¬ 
lapse” in the following, recognizing also that higher-order 
dynamics will play a role. 

To study stability, one is interested in the e —>■ 0 limit. 
In this case, for collapse to occur, spacetime quantities 
must continue to be large in this limit. Recalling that 
spacetime quantities are generally given as mode sums 
multiplied by powers of e, it would be necessary for these 
mode sums to diverge to see an indication of collapse. 
(We are supposing that jmax —t oo for this discussion.) 
For example, (j) = with (j)A'> given by ([^, so the 

only way for (j) to become large in the e —> 0 limit is 
for the sum ([^ to diverge. In this scenario it is possible 
to have perfectly well-defined TTF evolution, but with 
spacetime quantities ill-defined for any value of e. 

For exponential spectra, Aj ^ , sums such as 

always converge. But for power laws, Aj ~ (l-|-j) “, this 
is not the case. Indeed, at the origin, where the mode 
functions peak. 


ej(0) 


Wf + 3j + 2 
■ 


o{j). 


(47) 


So, for example, 

OO 

<^(t,0) = (A,(r)e—A + A,(r)e“A) e^.(o) 

j=o 

~e5](l+i)-“xO(j), (48) 

j 


^ In general relativity, collapse usually occurs once TTF breaks 
down [7|, while in Gauss-Bonnet gravity it can be averted m 
because of a radius gap for black hole formation. This holds 
despite both theories having identical TTF equations m- 


thus for a < 1, <() is UV-divergent. (We have here 
assumed that the phases do not cause a cancellation.) 
Other quantities, such as the metric variable A, are even 
more divergent. We therefore propose that the large-j 
asymptotic behavior determines whether black hole col¬ 
lapse can occur in the e —>■ 0 limit. (See also (30] for 
further discussion on this point.) 

Connecting to our study of QP solutions, the picture 
that emerges with regard to collapse is as follows. QP 
solutions that have asymptotically exponential tails will 
not collapse because they are equilibria and have well- 
behaved associated spacetime quantities. Initial data suf¬ 
ficiently close to a stable QP solution (with the same E 
and N) will also not collapse because the solution will 
simply oscillate around that QP solution, and its high-j 
tail will be close to the QP tail. Initial data that oscillate 
about a stable QP solution, but whose oscillations are 
quite large can collapse if the oscillation passes through 
a power law that causes the TTF description to break 
down. Finally, initial data that do not oscillate about 
QP solutions can attain a wider range of configurations, 
and, as we will confirm, tend to approach power laws and 
collapse (in AdS 4 ). 

In the following, we will examine several example so¬ 
lutions within TTF, both non-collapsing and collapsing. 
For the non-collapsing examples, our approach is to iden¬ 
tify the closest stable QP solution, and relate the ob¬ 
served dynamics to the linearized analysis. We find that 
the linearized eigenfrequencies {A„} (and combinations 
thereof) do a remarkable job of approximating the re¬ 
currence times, even for large perturbations. It should 
be kept in mind that, while the physically relevant limit 
takes Jmax —> c», all simulations are by necessity per¬ 
formed at finite jmax < oo; we will discuss the continuum 
limit below. 


A. Nearly-QP initial data 

We first study the nonlinear dynamics of initial data 
that closely approximate a stable QP solution. We show 
that the simulation closely matches the linearized anal¬ 
ysis, and we identify the origin of deviations from linear 
behavior. 

In anticipation of the following subsection, we dehne 
(a particular case of) two-mode initial data, 

= + (49) 

with the energy evenly divided between the two low¬ 
est modes. This data has temperature T = 3.75, and 
for later comparison with spacetime simulations we take 
E = 0.0162. There is therefore only one associated QP 
solution, with jV = 0 (see Fig. [^. (We neglect QP solu¬ 
tions that skip over modes.) Following m. we consider 
initial data that interpolates between the two-mode ini¬ 
tial data and the associated QP solution, 

Ej={l- X)Ef^ + (50) 
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FIG. 8. Evolution of mode j = 0 in complex plane for in¬ 
terpolated initial data. The full solution here is Ajij) = 
Note that here (and subsequently) 
we also fit for a A-dependent frequency shift satisfying /3j' = 
Po + j(/3i ~ Po)y which arises from nonlinear effects (it is 
quadratic in A). Had we not done so, there would be an addi¬ 
tional overall phase oscillation. This phase, however, has no 
influence on the evolution of the energy spectrum, and tends 
to 0 as A —>■ 0. 


where is the associated QP spectrum. For all A, this 
interpolation preserves E and N. 

We performed nonlinear evolutions of the TTF equa¬ 
tions ([^ for initial data (50) as the parameter A was 
varied between 0.05 and 0.30. Fig. shows the oscilla¬ 
tions of the lowest mode. As A is increased, the mode 
continues to oscillate about the QP solution, although 
with larger amplitude, as expected. We plot the evolu¬ 


tion of the energy of mode j = 5 as A is varied in Fig. 9a 


This shows that as the amplitude fluctuations increase 
considerably, the periodicity is not significantly changed. 
The discrete Fourier transform in Fig. |9b| shows that the 
oscillations are described by a discrete set of frequencies, 
as expected from the linearized analysis. 

Fig. [To] shows the peaks of the spectral energy density 
of 3?(Xj) for j < 20. For the most part, these peaks align 
closely with linear eigenfrequencies of the QP solution, 
but there are several extraneous peaks at low frequen¬ 
cies. These arise mostly in modes j = 0,1 and for larger 
A, which indicates they arise nonlinearly. This is con¬ 
firmed in Fig. El which shows that the new peaks grow 
nonlinear ly with A. In fact, the frequency of the first 
new peak is precisely the difference between the two low¬ 
est eigenfrequencies iXo = 0.0207 and iXi = 0.0396, so it 
is a nonlinear effect driven by a coupling between the two 
lowest eigenmodes. More generally, given the form (44) 
of the spectrum, 


iXji — Cl C2n O 1 

these lowest-frequency quadratically driven oscillations 
will arise at frequencies that are approximately C 2 = 
0.0158. At higher nonlinear order, additional low- 



FIG. 9. Energy content of mode j = 5 as a function of time 
(top), and its spectral density (bottom). 


frequencies will appear, at, e.g., C 2 — Ci = 0.0048. 

Notice also from Fig. [TOj that larger-j modes are influ¬ 
enced more strongly by the larger-n QP eigenmodes, as 
expected from Fig. Moreover, as the deviation from 
the QP solution increases (larger A) higher frequency QP 
eigenmodes are excited. 

This analysis shows that for nearly-QP initial data, the 
linearized analysis of the associated QP equilibrium so¬ 
lution does an excellent job of predicting the nonlinear 
dynamics, in particular the periodicities. Furthermore, 
as A is increased further, an additional low-frequency 
(« C 2 ) mode is nonlinearly excited by the linear oscil¬ 
lations. This mode, we will see, is most closely related to 
recurrences. 


B. Two-mode equal-energy initial data 

Setting A = 1 in the interpolated initial data of 
the previous subsection, we obtain the two-mode equal- 
energy initial data, which have received significant at¬ 
tention miiin]. As above, T = 3.75, so there is a single 
associated QP solution (that of the previous subsection). 

Fig. |12| shows the nonlinear evolution of the energy of 
the first six modes. The main recurrence time is closely 
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FIG. 10. Main oscillation frequencies present in the evolu¬ 
tion of nearly-QP initial data (for various A). Horizontal blue 
dotted lines represent the linear eigenfrequencies {| A„ |} of the 
QP solution. For each mode j, the three most dominant peaks 
of the spectral energy density are indicated by circles (largest 
peak) and crosses (secondary peaks). We also dropped sec¬ 
ondary peaks if they were smaller than 1% of the main peak. 
These plots were computed using a discrete Fourier transform 
on a regular grid of size 10'*. 


related to the periods of the QP eigenmodes, but is in 
fact slightly longer. Indeed, the dominant time scale of 
the 3 = 1 mode is approximately 450, and the first three 
linearized eigenmodes about the QP solution have peri¬ 
ods 303, 159 and 110. In precisely the manner described 
in the previous subsection for smaller A, nonlinear cou¬ 
plings between the eigenmodes drive oscillations at the 
new (slightly longer, as compared to the largest eigen- 
period, 303) characteristic time scale ~ 398 (for 

jmax = 100). Notice also that at the third recurrence 
(r « 1350) there is an even closer return to the initial 
configuration, and that this coincides with a third order 
nonlinear interaction time scale, 27r/(C'2 — Ci) « 1310. 

To compare with numerical simulations in AdS, we re¬ 
construct the spacetime fields from the TTF variables 
{Aj(r)} (with e = 1 by convention to keep time axes 
consistent). We plot the upper envelope of iP = 
at the origin, a: = 0, in Fig. 13 (11^ itself exhibits a fast¬ 
time oscillation that is not of interest.) This quantity is 
related to the Ricci scalar, and is frequently employed 
as an indicator of collapse (e.g. [illHlET]). Notice that 
{x = 0) can reach very large values in the course of evo¬ 
lution, but inherits the recurrences from the energy plot. 
Growth in n^(a: = 0) reflects direct turbulent cascades 
of energy to high-j modes, while decay reflects inverse 
cascades. The time scale of these recurrences—troughs 
at t = 450,850, peaks at t = 300, 700,1100—is consistent 



FIG. 11. Spectral energy density of 5R(xi), rescaled by A^. As 
A is varied, the heights of the two largest peaks are unchanged 
to leading order. Meanwhile, the rapid growth of the smallest 
peak (inset) with A indicates a nonlinear origin. Note also 
that Ai — Ao = 0.01881, which closely matches the position of 
the smallest peak. Thus, we conclude that the smallest peak 
arises from a quadratic coupling between modes n = 0,1. 
Note that we also observed even smaller peaks at frequencies 
Ai -I- Ao and A 2 — Ai. These plots were computed using a 
discrete Fourier transform on a regular time grid with step 
size 0.25, up to time 150,000. 


with the predicted period of 271/(72 ~ 398. 

It is now clear that previously unexplained recurrence 
times can be understood naturally as oscillations about 
QP equilibria, and they can be predicted without any 
time-integrations. (In the case of the two-mode data, 
the frequencjQ of recurrences emerges nonlinearly as 
the asymptotic separation C 2 between eigenfrequencies.) 
Such predictions are of particular relevance for their holo¬ 
graphic implications for field theories. 


C. Gaussian initial data, a — 4/10 

Initial data with a Gaussian distribution for the scalar 
field in position space have been closely scrutinized 
within the context of the AdS stability problem (see, 
e.g., [3 [13 HZ]). In particular, collapse was first stud¬ 
ied for a Gaussian with variance a = 1/16. As noted 
in |8], there is also a range of cr for which collapse is 
apparently averted. Armed with our new understanding 
of perturbations about QP solutions, we here analyze a 
non-collapsing Gaussian, and in the following subsection 
we study the collapsing case. 

The a = 0.4 Gaussian, which has T = 3.42, is in 
many ways similar to the two-mode initial data. There 
is a single associated QP solution, and the evolution 


In contrast to the frequency. Fig. |l3| shows that the amplitude of 
recurrences depends strongly upon jmax for the range we stud¬ 
ied. Since the oscillation is nonlinear, there is no obvious way to 
predict this amplitude. 
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FIG. 13. Upper envelope of Il^{x = 0) for two-mode initial 
data. Approach to full numerical relativity simulation is seen 
as j’niax is increased. The higher peaks arise because increasing 
imax allows the direct cascade to proceed to higher-j modes, 
which are more peaked about x = 0. This figure updates a 
similar figure in our previous work m, with a higher resolu¬ 
tion GR simulation, and also larger jmax TTF simulations. 


be determined from the linear spectrum to be 12500 for 
jmax = 100, in agreement with Fig. 
also matches the recurrences in Fig. 


This time scale 


FIG. 12. Energy evolution of first 6 modes for two-mode 
equal-energy data with T = 3.75. The evolution was com¬ 
puted using the TTF equations with jmax = 100. 


is characterized by a series of direct and inverse cas¬ 
cades. Throughout the evolution, the energy spectrum 
(see Fig. remains roughly exponential—as opposed 
to power-law—corresponding to non-collapse. 

Observed oscillation periods can be predicted by an¬ 
alyzing the associated QP solution in the same way as 
for the two-mode data. As a representative example, we 
monitor the behavior of a high-frequency (j = 58) mode 
in Fig. 15 The high-frequency oscillation of 3?(x58) oc¬ 
curs with period 359, which to this accuracy matches 
exactly the period of one of the A„. Similar agreement 
with the linearized frequencies can be seen for the other 
modes. 

The Gaussian data, however, differs from the two¬ 
mode initial data in that, initially, it more strongly ex¬ 
cites high-j modes. In turn, this causes increased ex¬ 
citation of high-n QP eigenmodes. This is reflected in 
complicated linear “beating” and nonlinear “driving” dy¬ 
namics between excited modes seen in Fig. |15[ Here, 
the slow envelope modulation arises as the difference in 
frequencies between subsequent QP eigenmodes—with a 
corresponding period l-KjCi for this QP solution. The 
amplitude of the beating is predicted by the linear anal¬ 
ysis to be ^ 10% of the measured amplitude and we ex¬ 
pect that nonlinear driving accounts for the remainder. 
Again, the characteristic time scale is 27r/C2, which can 


D. Gaussian initial data, a = 1/16 


In contrast to all previous examples, the cr = 1/16 
Gaussian is seen to collapse in numerical simulations [71 
IS]. The temperature T = 13.1 suggests that there could 
in principle be several associated QP equilibria, but nev¬ 
ertheless the data do not display any oscillations, indi¬ 
cating that they are far from these equilibria. 

Consistent with previous full GR simulations |S3], the 
enercy spectrum of this data in TTF approaches a power 
lawnijj ~ (j -I- 1)““ as it evolves in time (see Fig. 161. 
Extrapolating to jmax oo, such a spectrum would lead 
to diverging spacetime helds (such as H^ at the origin), 
indicating the break down of TTF as a valid description. 
At this point, higher order dynamics have been seen to 
lead to collapse [ZlIH]. Despite the failure of TTF to pro¬ 
vide a valid description past the power law, the TTF so¬ 
lution (for finite jmax) is perfectly well-defined for longer 
times (beyond those shown in Fig. 161. 


A recent publication m examined the TTF evolution 
of two-mode data in AdSs, which displays a similar evolu¬ 
tion to a power-law spectrum as seen here. It was argued 


® Speculation [30] that power laws do not arise in gravity (in an 
analogy to a self-interacting scalar field) are based on scaling as¬ 
sumptions for the tS-coefficients. In fact, the tS-coefficients grow 
with increasing mode number in gravity (see footnote |10[ |, while 
they decay for the scalar m, so the coupling to high modes is 
much stronger in gravity. This arises because of the spacetime 
derivatives present in the gravitational interaction. 
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FIG. 14. Evolution of the energy spectrum for a = 0.4 Gaus¬ 
sian initial data. We show several times during the first direct 
and inverse cascades (top), and a much later time during the 
second inverse cascade (bottom). Spectra are all roughly ex¬ 
ponential. (jmax = 200) 

that in the limit jmax —>■ oo, TTF itself breaks down af¬ 
ter the power law is reached, with the time derivatives 
of phases of the mode amplitudes oc log(T — t*). The 
truncated equations nevertheless have a well-defined so¬ 
lution beyond t», which was described as an unphysical 
“afterlife.” The authors of m emphasized that even at 
finite jmax a highly oscillatory behavior led to numerical 
difficulties beyond r*. We note that we did not encounter 
any numerical difficulties in our simulations of ct = 1/16 
Gaussian data beyond this tim^^ (see App. 0 - 

We present in Fig. the evolution of the first deriva¬ 
tive of the phase of mode 80 for different values of /max- 


The closed-form expressions of App. ^ give rise to analytic ex- 
pressions for the asymptotic scaling of the tS-coefhcients. Pre¬ 
cisely, in AdS 4 , ~ In j. For comparison, Ref. |31| 

reports the corresponding expression in AdSs to be j^i^. Be¬ 
cause of the Inj- factor, the arguments of m do not apply in 
AdS 4 ; the phases (in their notation) (Bn) ~ nlnn can have 
Bi -h Bn — Bj — Bk oo even for resonant quartets, so the ansatz 
taken for high modes is invalidated (e.g., B 2 i-\-Bo — 2Bi ~ 2 In2). 
Since it would be natural for a logarithmic factor to arise in AdSs 
as well, it would be useful to analytically compute in this case 
the asymptotic form of the tS-coefficients. 



xlO”® 


T 



0 10000 20000 30000 40000 

T 

FIG. 15. The evolution of 5R(x58), and the upper envelope 
of n^(a: = 0), for a = 0.4 Gaussian initial data. The high- 
frequency oscillation of 1 R(X 58 ) is very well predicted by the 
linear analysis of the associated QP solution. The lower- 
frequency modulation is a combination of beating of linear 
modes and nonlinear driving, and it corresponds to the recur¬ 
rences seen in 11^ (a; = 0). (jmax = 200) 


While the curves agree at early times, we were not able to 
conclude whether they approach a limit near the collapse 
time T* « 1500 as jmax —^ oo. Nevertheless, we do not 
observe the logarithmic blowup in the derivatives of the 
phases reported in Fig. 5 of [3T] for the case of AdSs; in 
the present case of AdS 4 the behavior seems less extreme. 


V. CONCLUSIONS 

In this paper, we have analyzed perturbations of AdS 4 
within the TTF formalism. We identified a collection 
of two-parameter families of QP equilibrium solutions 
to the TTF equations, and we established their linear 
stability (in the sense of Lyapunov). For each QP solu¬ 
tion, this analysis gave rise to a new spectrum of eigen- 
modes {e„}—which are collective oscillations of the AdS 4 
normal modes {ej} about the QP solution—along with 
their own oscillation frequencies {A„}. We also showed, 
through several examples, that the linear analysis often 
remains valid well into the nonlinear regime, and more¬ 
over, the leading nonlinear effect is generally to introduce 
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FIG. 16. Energy spectra at various times for a = 1/16 Gaus¬ 
sian initial data. This figure shows the approach to power- 
law during the time evolution of the TTF equations. We 
used /max = 400. Gompare to the exponential spectra for the 
(7 = 0.4 Gaussian, illustrated in Fig. |14| 



FIG. 17. Absolute value of the derivative of the phase of mode 
80, for an initial Gaussian o = A data, and different values 
of /max. We denote Aj = 


yet the QP solution [in conjunction with conservation of 
{E, N, H)] constrains the available phase space and acts 
as an island of stability. The spectra of perturbations 
about QP solutions are non-resonant themselves (except 
asymptotically for large n), indicating that thermaliza- 
tion should not be expected within a stable island. 

We point out two facts that play a role in extending 
this work to systems beyond AdS 4 . In AdS, each QP 
family is parametrized by the two conserved quantities 
E and N. In contrast, for other confining systems, such 
as a flat spherical cavity [53], additional resonances are 
present, and N is not conserved. As a result, in that case 
we expect only one-parameter families of stable equilib¬ 
ria. Meanwhile, as the dimensionality of the system is 
increased, couplings to high-j modes become stronger 
(5-coefficients become larger), which then drive stronger 
turbulent cascades. These effects compete in deciding 
collapse versus non-collapse. 
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new frequencies that are combinations of the {A„}. 

A key takeaway message is that for initial data that do 
not collapse as e —>■ 0 in AdS 4 (or, at least, do not do so 
immediately), recurrences are simply oscillations about 
stable QP equilibria. The relevant frequencies arise from 
the {A„}. With our stability analysis, we now have a 
method of predicting recurrence times without any need 
for time integrations. 

Initial data that do collapse as e —>■ 0 (in the sense that 
the TTF description breaks down) are not sufficiently 
close to any QP solution. We observed, in agreement 
with previous fully nonlinear general relativistic simu¬ 
lations [24], that these data tend to approach power- 
law energy spectra. The presence of stable QP equi¬ 
libria thus reconciles the apparent tension between the 
fully commensurate frequency spectrum of AdS and non- 
thermalizing initial data. Indeed, fully commensurate 
frequency spectra would be expected to thermalize as 
a result of the KAM theory and Arnold diffusion [55] . 


Appendix A: Closed form for the iS-coefHcients 

Reference El provided new simplified formulas for the 
5-coefficients described by integrals of products of the 
mode functions 0 . Notice that our conventions for the 
5-coefficients differ from those in El by a factor 4, that 
is to get our coefficients Sijki one has to multiply the 
expressions given in [14] by 4. Note also that in this 
section, for commodity, we rewrite Siju = Sjki- 

Here we shall give new closed-form formulas found for 
the tensors of which allowed us to compute the 5- 
coefhcients up to /max = 400 in a short time. 

Recall the expressions for Sijki given in [15] : 

Sun = 2LofXuii + 6Yiui + SojfWiiu + 8ujfWi*ui 

— 4ujf{Aii + ujfVii), (A-1) 
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and, for i I, 


^lili — — 2( 


'^l -^i 


){ujfXini — ^iXiiii) 


+ - Xuii) 

ujf - ujf ujj - w2 

+ 2{Yuu + Yuu) + AuJfoJf{Wu^^ + Wuii) + 

+ ujfW^^ii — tJf{Aii + ujfVii). (A2) 

and finally, i ^ I and i ^ 

^ijkl ~ 

( I “ 1 “ “ 1 “ '){j-^i^j^k^lijk ^l^iljk') 

bJi + UJj UJi — UJk OJj — ^k 

( \ “ 1 “ j^k^l^ijkl ^iA^jikl) 

LUl -j- UJj UJi — UJf^ UJj — UJi^ 

( I “t“ ^{y^i^k^lXjilzl 

UJi -j- UJj UJi — UJk UJj — UJ}^ 

(tT I 71 7^ 77 ^^'ji^i^j^lXkijl ^kYikjl) • 


Since the first argument of the hypergeometric function is 
a negative integer, the hypergeometric function appear¬ 
ing in Cj is in fact a polynomial function of degree j. 
Using then the expansion for the hypergeometric func¬ 
tion, 


3 \ ,^fc(3 + j)fe r 


-J,3 + j;-;xU^ (-1) 




' a)k 


(AlO) 


with {a)k the rising Pochhammer symbol, the following 
identity holds as proven below: 


ej{x) = 




\/<j + 1)(J + 2) sin(a;) 


(All) 


with 


fm{x) = (m + 1) sin(2TOx) + msin[2(m -|- l)a;]. (A12) 

To establish this result, first notice that in d = 3, the 
mode functions satisfy the differential equation, 


I j UJi k ^^3 k 


(A3) + 2 [tan(a:) + cot(a;)] e'(x) + (2i -I- 3)^ei(a;) = 0. 


The quantities that appear in these coefficients are de¬ 
fined by integrals of the mode functions (recall that we 
work here in d = 3 spatial dimensions): 

f ^ sin^ X 

Xijki = / ei(x)ej{x)ek{x)ei{x) -da;, (A4) 

do cos a; 


E'fc(a;)e;(a;)^^: 

^ ^ cosa; 


(A13) 


Next, using (All) and ([^, it is immediate to check that 
at a; = 0 the two expressions and their derivatives have 
the same limits, 


(A4) 

hm ei(x) = 4U 

lim e'Ax) = 0, 

(A5) 

lim e"(a;) = -- 

fc—vO 


(* -f l)(f + 2) 


(A14) 


Wail = / da;ei(x)^ sina;cosx / dj/efc(?/)^ sin j/cosy, 
Jo Jo 

(A6) 


^iiii — / da;e'(x)^sina;cosx / dy efc(y)^ siny cosy, 
do do 


(A7) 


Vij = / da;ei(a;)ej(a;) sina:cosx, (A8) 


Aij = / dx e((x)e'(x) sinx cos X. (A9) 
/o 


To simplify these expressions, we used the form of the 
mode functions ej{x). We know that ([^ 


ej(x) = 4:\j cos^ X 2 T’i ( — j, 3 + j; sin^ x 


It is then straightforward to check that both (All) 
and ( |A10 ) satisfy the differential equation (A13) on 
(O, At this point, while one might be tempted 
to conclud e bo th expressions are the same, we 
note that (A13) is singular at x = 0 and x = 


We thus proceed as follows: let us first denote 
ej(x) = cos(x) u, (sin ^ x). Next, with the expanded 
form of both ( |A10 ) and ( All[ ) [using sin(2mx) = 

show that Uj(t) is in both cases a polynomial function 
of t; the remaining task is to show both polynomials are 
the same. To check this fact, denote by \Q{t), T{t)} 
the polynomial equal to u(t) from expressions (AlO) 
and dAlTt , respectively. We ca n the n substitute each 
(multiplied by cosx) into Eq. (A13|. The resulting 
equation for both Q and T in (0,1) is the same simple 
differential equation. 


+ (^^(1 — t) — 3 -I- = 0, 


(A15) 


with / standing for either Q or T. Now, since Q and 
T are polynomials, they verify this equation everywhere. 
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Moreover, this equation gives rise to the following order-2 
relation on the coefficients of Q and T, denoting Q{X) = 
and T{X) =J2 which is 

k k 


2{k + l)(2k + 3)uk+i + (w| — 3 — — 6k)uk 

+ (4fc(fc — 1) -I- 1 — ujj)uk-i = 0. (A16) 


Consequently, according to (A16), {qk)k and (tk)k are 
both uniquely determined by the same relation, and by 
their first two values. It is thus sufficient to check that 
qo = to and qi = ti, which is given by (A14) [since e(0) = 
u{0) and e"(0) = —u{0) -f 2u'(0 )1. We have thus proven 
that Q = T, that is that (All) is a valid expression for 
ej- 

Let us also stress that these calculations, done in AdS 4 , 
are not straightforwardly extended to other dimensions. 
Some inspection and analysis of the calculations in dif¬ 
ferent dimensions, however, points towards a similar sim¬ 
plification of eigenmodes in odd spatial dimensions d, 
though we have not exhaustively studied this question. 

We then have, for the derivative, 


e'Jx) = 


3 + 2j 


cos(a;) 


+ 1 ) 0 ' + 2 ) sin^ix) 


9 j+iix), (A17) 


with 

gm{x) = —{m + 1) su\{2mx) -\- msin[ 2 (m -|- l)a;]. (A18) 

Since the indefinite integrals appearing in W 
and W* are easy to compute, one can now re¬ 
duce the problem to computing many integrals 
of the type da: a;'*'cos“ a; sin^ a; x F{2mx) and 
/p^ da; a:'*'cos“ a: sin^ a: x F{{2m -\- l)a:), with F the co¬ 
sine or the sine, m an integer, 7 S { 0 , 1 }, and a and /3 
integers greater or equal to — 1 . 


We give here some conventions and the few delicate 
integrals that one has to compute in order to get the 
relevant coefficients. We will denote 5^ the Kronecker 
delta function, and Sign(m) the function taking the value 
1 on N*, the value —1 on Zl, and Sign(O) = 0. Some 
of these integrals were found thanks to several formulas 
found in [51] . 

We also denote V' the polygamma function ^l){x) = 
where L is the Euler function. We denote p = m + n 
and k = m — n: 


Vn S Z, 


, sin ((2n -I- l)a:) 7 r,„. , , _ ,, 

da;-—--= -(Sign n) -b <5 n)), 

sm (x) 2 


Vm € N, 


dxx °"((^’"+i)") = ^(-ir 

sm(a;) 4^ ^ 


V(to, n) G 


' m + 2 lA + m 


, cos (a;) . N . /r, X 

da: ^sm (2mx) sm(2nx) = 
sm (x) 


if p is even 

-V’(f) + ^(|) + i-^] if pis odd 


We are then able to find closed form expressions for every 
quantity we need. Nevertheless, these expressions appear 
to be too long for the X and Y tensor to be written clearly 
on one page, and are therefore not given here. We shall 
now give the expressions found for A, V, W and W*. 


Due to the symmetric property of A and V it is sufficient to restrict to i > j. For the case where i + j is an even 
integer. 


Aij — 


9 /(■ + + + + 

27r^(? + !)(? + 2)(j + l)(j + 2) V 

2(i + l)(j + 1) {-2{i - - 2i{i + 8 )j + i{2i{i + 4) - 3) + 2j^ - 3j - 


* + j + 3 


2 

13)' 




(A19) 


27r \/(f + 1)(* + 2)(j + l)(j + 2) 


(2z + 3)(2j + 3) 


' A + j + lx + 

V'(-^-) - 'iPi —^—) 


(2f-b3)2 i{i + 2) 


(i -b l)(i -b 3) 6 — 8i{i - 

—i+j — 1 i+ j + l 


1 ) 


■4(3i-b4) 


-b 


i+j + 3 i-j-1 

(A20) 
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For the case where i + j is an odd integer, and i > j, 


Aij — 


(2i + 3)(2j+3) 


27r\/(* + 1)(* + 2)(j + l){j + 2) 


(2.(. + 3) + 2,U + 3) + - V'(^)) - 4.; - ' 


* - J 


7z(i + 2)+6 (i + l)(i + 3) 


i+j+ 2 


i + j + 4 


- 8i 


(A21) 


= 


47r(i - j)A/(i + l)(i + 2){j + l)(j + 2)(i + j)(i + j + 2)(i+j + 4) 
8(i((i - 21)i - 93) - 89)f + 4(i(i(i(6i + 37) + 19) - 102) - 108)j 


(-8(5i + 7)j^ - 4(i(14i + 85) + 93)^^+ 


-2(2^ + 3)(2j + 3)(^ - j)(* + J)i^ + j + 2){t+j+ 4) 




+ 16i(i + l)(2i(i + 5) + 9) . (A22) 


Last, we have the following general values for W and W*: 

{-2m'^{2l + 3) + Am{P - 2) + 21{ZI + 5) + 3) 




167r(m + l)(m + 2)(l + l)(l + 2) (21 + 3) 

4^mTlKm +2) +!) + >/-("'+ j) + 21»(2)) 




(2m + 3)2(2Z(Z + 3) + 5) 


167r(TO + 1 )(to + 2){l + 1)(; + 2){2l + 3) 


Sign(m — 1) 


1 


[8m^(l + 1)(2;(Z + 4) + 7) + 8m^{l + 1)(Z(14/ + 55) + 48) 


167r(m + l)2(r7T, + 2)^(Z + 1)(Z + 2)(2l + 3) 

+to2(4Z(;(73Z + 355) + 527) + 979) + + 24) + 602) + 1113) + 2{1{1{1AI + 351) + 515) + 237)] , (A23) 


W 

r r n 


(2m + 3)^ 


mmnn 


167r(m + l)(m + 2)(n + l)(n + 2)(2n + 3) 

+4m(n(n(6n(n + 3) — 13) — 54) — 18) + 3(2n(n(6n^ + 28n + 41) + 21) + 9)] 5m- 


[4m^(2n + 3) — 8m^((n — 3)n — 7) — 2m^(2n(6n(n + 6) + 41) — 3) 


(m + 2)(2m + 3)^(n + 1)(—m + n + 1)^ 


5m—n—l 4 “ 


(m + l)(2m + 3)^(n + 2)(m — n + 1)" 


47r(m + l)(n + 2)(2n + 3) 7r(m + 2)(n + l)(8n + 12) 

(2m + 3)^ (4TO^(2n(n + 3) + 5) + 12m(2n(n + 3) + 5) — 2n(n + 3)(8n(n + 3) + 27) — 37) 
167r(m + l)(m + 2){n + l)(n + 2)(2n + 3) 

(2m + 3) 


5m—n+l 
’Sign(m — n) 


^[l6m^(n + l)(2n(n + 4) + 7)+8m4(n + l)(2n(17n + 67) + 117) 
167r(m + l)^(m + 2)^(n + l)(n + 2)(2n + 3) 

+m3(8n(n(n(93 - 4n) + 517) + 797) + 3018) + m2(4n(n(n(187 - 36n) + 1458) + 2408) + 4701) 

+m(4n(n(n(31 - 52n) + 936) + 1736) + 3545) + 2(n(n(423 - 2n(24n + 31)) + 947) + 513)]. (A24) 


Appendix B: Methods for obtaining quasi-periodic 
solntions 


r 


find numerical QP solutions of (141 it is necessary to 


choose an appropriate “seed” for the algorithm. 


Here we describe the two approaches we took to gen¬ 


erate the families of QP solutions in Sec. IIB 


1. Direct solution using Newton-Raphson method 

Given an appropriate starting point, the Newton- 
Raphson method provides successively better approxi¬ 
mate solutions to a set of coupled equations. Thus, to 


Although we parametrized QP solutions by E and N 
in the main text, it is more appropriate here to choose 
parameters from among /Sq, Pi and {aj}. For example, 
to find jr = 0 solutions we fix ao = 1 (an arbitrary choice 
because of the scaling symmetry) and ai ^ ao. Eqs. (14) 
may be solved to eliminate Po and Pi, leaving jmax — 1 
equations and jmax ~ 1 unknowns. For the remaining 
variables, we choose an exponential energy spectrum as 
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a seed, 


2 J+ 3 ’ 


(Bl) 


with /i = log [3/(5ai)]. For sufficiently small a\ the 
Newton-Raphson method gives back a solution with 
nearly exponential energy spectrum, as seen in Fig. 

As ai is increased, the QP energy spectra deform away 
from exponentials and it becomes increasingly difficult 
to find solutions with the exponential ansatz ( |B1[ ). In 
fact, in [T2] we could not find solutions with oi > 0.42. 
Slightly better results can be obtained by taking (/3i — 
/lo)//3o as a parameter {Pi » Pq approaches the single¬ 
mode solution), however, this also breaks down for large 
T. To fully uncover the jV = 0 family we require the 
technique of the following subsection. 

Solutions within jr > 0 families can be obtained in a 
similar manner, now fixing aj^ = 1 and ^ ocj^. 

To pick a seed for the Newton-Raphson algorithm, we 
solve the first several QP equations (14) perturbatively 
in . 


2. Perturbation from known solution 

Now suppose is a known numeri¬ 

cal QP solution. Sec. |IIIB 2| shows that there is in gen¬ 
eral a 2-parameter family of perturbations to nearby QP 
solutions, so by following such perturbations new QP 
solutions—otherwise not readily obtainable through the 
Newton-Raphson method—can be constructed. 

Since one of the parameters is, as usual, an overall 
scale, there is only one nontrivial parameter. It is, there¬ 
fore, convenient to fix N and vary if by a small amount 
6E. Following Sec. |IIIB^ we numerically solve the linear 
system of equations, 


0 = 2ijjj [aj{9i + ujj92) + PjUj] (B2) 

T ^ ^ ^klm 4“ 4“ ^ 

klm 

0 = 8^a;jajUj, (B3) 

3 

SE = 8 a;| aj Uj , (B4) 

3 

for the variables (0i, 02, {wj})- We then update the QP 
solution. 


—>■ q; j 4-Mj, (B5) 

Pj —y Pj 4- 01 4- ujj92- (B6) 


Newton-Raphson method can be used periodically to en¬ 
sure the deviation from actual QP solutions does not be¬ 
come too large.) In this manner, we obtained the full 
QP families illustrated in Fig. [3| These families termi¬ 
nate when solutions to (B2)-(B4) no longer exist (i.e., 
when the associated matrix has vanishing determinant). 


Appendix C: Minimization of H and linear stability 
of QP solutions 


We shall here show the relation between the minimiza¬ 
tion of E[ for a QP solution and its linear stability. We 
know from (22) that has a critical point at a QP so¬ 


lution. Here we compute the second order change in iJ. 

Let us take a generic second order perturbation of a 
QP solution, that does not perturb E and N, 


Q! '} 




( 2 ) 


(Cl) 


ik') 

where A\ is the order k perturbation of aj. Recall the 


expression 0 of H, 


1 


^ ^ ^ ^klm.^^3 

jklm 


A.AkAiAr 




Inserting © into this equation, one finds 


E 




l(l)|2 


4-^ ^ ajakaiA^^^ + ajakamAj^^ + ajaiamA 

j,k,l,7n 


( 2 ) 


+akaiamA^‘^^ + aiamAf'’ 
EoikamAj ^ + akOiiAj 
+ajaiAl}'’A^^'> 




Let us concentrate on the part where only Aj ' appears. 


Using the QP TTF equation (14), as well as the rela¬ 


tions (15) and (16) on the S coefficients, one can reduce 


the expression of this part to 

X! “4 “ X! ^3^3Pj )■ 

fc 3 3 

(C2) 

Now, since E and N are conserved at both linear and 
quadratic level, we have, for the quadratic level. 


E 

3 

E 




aj{A 




( 2 ) 

3 

( 2 ) 


■Af) + \A 


( 1)|2 
3 I 




= 0 , 


= 0 , 


The new QP solution has particle number A, energy E + 
SE, and therefore the temperature has changed by ST = 
SE/N. 

With the updated QP solution, the procedure may 
be iterated repeatedly to obtain finite-sized AT. (The 


which can be rewritten as 

y{uj) G s.t. Uj = uq + j{ui — uo), 


E‘ 


^3 ^3 


,{Af^+Af>) + \A 


(2)^ 


1‘‘V 


= 0 . 
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Using this identity in (C21, one can rewrite the full second 
order variation of the Hamiltonian as a function of the 
linear perturbation only, 




k j 



jklm 


^ + OikOtlA^^ ^ ^ 

+ajaiA^^'^ 

Now, with ( [Tsl , one can reduce this last expression to 

- = f E I'+E 

3 3 

-^Ckal^ujjaj\Af\^ + ^ S^i^auaiAf^A^'^ 

k j jklm 

+ ^ ^ ^jklm0^l(^m[Aj Af, + Aj Af, ). 

jklm 


Let us now rewrite the linear perturbation A^j'^ in terms 
of real and imaginary part, 

SAj = Rj ilj- 

If we denote by X the column vector 
(i?0, • ■ • ) Rjmaa: ) -^0) • • • J )) ^>^^1 M thC ITiatriX 

such that we have —S^H = X'^MX, them M is of the 

, where A and B' are both square 

matrices of size jmax + 1, and we have the following 
expressions for their coefficients: 

2 _ ^ 

Im 

k kl 

Ri,j ~ ~ 2 E! 

Im 

- E + E ^fkiAkai. (C4) 

k kl 


simple form 


A' 0 
0 B' 


Note that these matrix elements are quite similar to the 
matrix elements of the matrix A whose elements can be 
deduced from ([^ and p9). 


Indeed, writing A in the form 


0 -C 
D 0 


one can, with 


the same type of calculations, prove the following simple 
identities: 


A^j = ujiD^j - 2aiaj{ujjCi - ujjCj), (C5) 

BX=oj,a,j. (C6) 


In (C5), since we are interested in the sign of X^A'X to 


characterize stability, the right antisymmetric part will 
play no role and we can ignore it. Let us also recall 
that we are interested in the sign of with X 

satisfying the linear conservation of E and N, that is, if 
^ {RQi ■ ■ ■ ) AOj ■ ■ ■ ^ Ijmax) ) 


y{uj) S S.t Uj = Uq + j{ui — Uq), 

'^j^ajUjjUjRj = 0 , 

j 


which is equivalent to saying that (Rj) is orthogonal to 
the vectors Xi = {ajUj) and X 2 = {ajuj) in the Eu¬ 
clidean Rtmax+i space. We will thus place ourselves in the 
two spaces E = {xi,X 2 )^ for A' and D, and E = ]RJ">ax+i 
for B' and C. 

We note that in order to get the announced result, one 
has to assume that C and D are both diagonalizable. We 
have seen numerically this is the case, but we have not 
rigorously proven this. 

Let us now assume that El has a local minimum at the 
QP solution (oj). Then that means that A' and B' are 
negative, 

VX e E, X’^A'X < 0, (C7) 

VF e E, Y^B'Y < 0 . (C8) 


/wo 0 0\ 

Denoting T = j 0 ... 0 , we have A' = TD. 

\ 0 0 w, / 

y 7max y 

Taking any eigenvector X of D with eigenvalue A (since 
D is diagonalizable), we have X'^AX = X'^ujiX'^ < 0, 

i 

which means that Sp{D) C K_. Using the exact same 
trick one can show that. 


TD is a negative symmetric matrix Sp{D) C ffi_. 

(C9) 

But, by computing the expression for the D coeffi¬ 
cients, it is immediate that DT~^ is also symmetric. 
Now, since we know that A = TD is a negative ma¬ 
trix, we deduce that DT~^ is also negative. Since TC is 
also negative, and since the non-zero eigenvalues of the 
product of two negative matrices are positive, the real 
eigenvalues of DC are all positive. Notice that since DC 
and CD have the same characteristic polynomial, this is 
also the case for CD. 

Let us recall that we argued that our system (28)- 


(29) is stable if and only if the eigenvalues of A are all 
pure imaginary. This means, since by deriving (28)-(29) 
again one can decouple the system of equations, that the 
eigenvalues of A^ are all real and negative. But since 
A2 ^ /-CD 

, we know that the spectrum of 


is going to be in IR_ if H has a minimum at the QP 
solution. So we know that if E[ has a minimum at a QP 
solution then this solution is linearly stable. 

Notice that this reasoning also holds if El has a maxi¬ 
mum at a QP solution; in that case the solution will have 
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as in EH, the implicit Runge-Kutta scheme of order 6. 
As an illustration, we present here two tests of the va¬ 
lidity of the implicit scheme adopted and our strategy to 
ensure no relevant short-time-scale is discarded. We held 
fixed the double-precision employed and varied the pre¬ 
cision of our adaptive step size method by 16p/10 digits, 
p = 1,..., 10. Fig. [T8| illustrates the change in conserved 
energy E vs integration time for different values of p. As 
is evident in the figure, the error quickly converges to a 
limit function, which is already reached for p = 4. (The 
remaining error is due to the double precision numbers.) 
We note that the results presented through the paper 
have been obtained with p = 5 with the implicit method. 


FIG. 18. The error in the total energy E for a Gaussian initial 
data of variance o = 0.4, using different values of p. We used 
imax = 100 for these calculations 


only unstable modes (we however never observed such a 
solution). 


Appendix D: Numerical integration method 


The integration of the TTF equations ^ requires spe¬ 
cial care as, depending on the values of the coefficients 
Aj , they can become stiff. Stiff equations can be handled 
with explicit methods—where the timestep must be small 
enough to ensure stability—or implicit methods—where 
stability issues can be more easily avoided but care must 
be exercised so as not to “discard” relevant short-time- 
scale physics by adopting too large a timestep. We have 
implemented both explicit and implicit methods as well 
as performed self-convergence in our analysis to ensure 
the correctness of the obtained results. 

In particular, we have employed the explicit (predictor- 
corrector) Adams method as well as backward differen¬ 
tiation formulas (both with adaptive timestepping) and, 


To further illustrate that no relevant short-time-scale 
physics was accidentally discarded by the use of an im¬ 


plicit integration scheme, we show in Fig. 19 the evolu¬ 
tion of a representative mode j = 50 mode for two rather 
distinct values of p. The difference between these two 
xlO”® 
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FIG. 19. The evolution of the real part of the 50*’’ mode, 
for the extremal values of p we took, 1 and 9. The difference 
between the two curves is of order 10”*”. We used jmax = 100 
for these calculations 
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